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The present generation of general circulation models (GCM) use parameterized cumulus schemes and 
run at hydrostatic grid resolutions. To improve the representation of cloud-scale moist processes and 
land— atmosphere interactions, a global. Multi-scale Modeling Framework (MMF) coupled to the Land 
Information System (LIS) has been developed at NASA-Goddard Space Flight Center. The MMF— LIS has 
three components, a finite-volume (fv) GCM (Goddard Earth Observing System Ver. 4, GEOS-4), a 2D 
cloud-resolving model (Goddard Cumulus Ensemble, GCE), and the LIS, representing the large-scale 
atmospheric circulation, cloud processes, and land surface processes, respectively. The non-hydrostatic 
GCE model replaces the single-column cumulus parameterization of fvGCM. The model grid is 
composed of an array of fvGCM gridcells each with a series of embedded GCE models. A horizontal 
coupling strategy, GCE <-* fvGCM <-> Coupler <-► LIS, offered significant computational efficiency, with the 
scalability and I/O capabilities of LIS permitting land— atmosphere interactions at cloud-scale. Global 
simulations of 2007—2008 and comparisons to observations and reanalysis products were conducted. 
Using two different versions of the same land surface model but the same initial conditions, divergence 
in regional, synoptic-scale surface pressure patterns emerged within two weeks. The sensitivity of large- 
scale circulations to land surface model physics revealed significant functional value to using a scalable, 
multi-model land surface modeling system in global weather and climate prediction. 

Published by Elsevier Ltd. 


1. Introduction 

The land and atmosphere form a highly coupled system. Surface 
heat and momentum fluxes are linked to the surface net radiation 
flux, the vegetation state, and the profiles of temperature and water 
from below the surface up through the atmospheric boundary 
layer. The fluxes of heat, momentum, and moisture across the land/ 
atmosphere interface are influenced by the heterogeneous char- 
acter of the land surface layer and vary on spatial scales ranging 
from meters to thousands of kilometers. Linking the water and 
energy cycles is precipitation. Feedbacks between the heteroge- 
neous land surface and the boundary layer affect the development 
of clouds and precipitation (review in Pielke, 2001). The vertical 
distribution of latent heat released through the formation of clouds 
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and precipitation modulates the large-scale atmospheric dynamics 
of the low and mid-latitudes, affecting the distribution, intensity, 
and longevity of waves, jets, and fronts, and thus to future 
precipitation patterns. Coupling a general circulation model (GCM) 
to a land surface model (LSM) allows for two-way interaction of 
atmospheric moist processes with the land surface. By coupling 
a GCM to a multi-model Land Information System (LIS) rather than 
to a single LSM, significant additional physical and functional 
flexibility is achieved (Kumar et al„ 2006; Peters-Lidard et al., 
2007). This paper describes the NASA-Goddard finite-volume 
Multi-scale Modeling Framework— Land Information System 
(MMF— LIS), a global model framework capable of explicitly 
resolving cumulus convection and simulating cloud-scale land/ 
atmosphere interactions. The MMF— LIS integrates an atmospheric 
GCM with a 2D cloud-resolving model (CRM) for explicit simulation 
of cumulus clouds and couples the LIS to the GCM. We describe the 
development and operation of the current Goddard MMF— LIS, 
focusing on the model coupling and its initial testing, particularly 
with respect to surface variables. This paper can be viewed as 
a third companion to two previous papers on LIS, the first 
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description of LIS by Kumar et al. (2006) and the role of LIS 
in coupled mesoscale modeling by Kumar et al. (2008). The 
MMF— LIS enhances our ability to investigate the integrated impact 
of small-scale cloud microphysics and soil and vegetation states 
on regional to global-scale circulations, cloud patterns, and 
precipitation. 

2. Background: global multi-scale modeling 

The current generation of GCMs used in operational global 
weather and short-term climate forecasting by the National Centers 
for Environmental Prediction (NCEP), the European Centre for 
Medium-Range Weather Forecasts (ECMWF), and the NASA Global 
Modeling and Assimilation Office have fully interactive land/ 
atmosphere coupling using single LSMs, respectively, Noah (Ek 
et al., 2003), Hydrology-Tiled ECMWF Scheme for Surface 
Exchanges over Land (H-TESSEL, Balsamo et al., 2009), and Catch- 
ment (Koster et al., 2000). Although the NCEP Global Forecast 
System (GFS) uses the NASA LIS for land data assimilation, only 
Noah is fully and directly coupled to the GFS atmospheric model 
component (Saha et al., 2010). These operational LSMs use tiles or 
catchment sub-divisions to improve the representation of the land 
surface heterogeneity within GCM gridcells. However, the surface 
fluxes generated are spatially averaged so that the atmospheric 
component can use a parameterized cumulus scheme to determine 
gridcell clouds and precipitation. Model comparison projects in the 
Global Energy and Water Cycle Experiment (GEWEX) have shown 
that simulations of various types of clouds and cloud systems from 
different geographic locations by cloud-resolving models (CRM) 
agree with observations better than those from cumulus parame- 
terizations used by the current generation of GCMs ( Moncrieff et al., 
1997; Randall et al„ 2003b). The lumping of land/atmospheric 
interactions and the use of cumulus parameterizations for cloud- 
scale moist processes are sources of significant uncertainty in 
predictions at larger scales (Zhang et al., 2005; Pauluis and Garner, 
2006; Shutts and Palmer, 2007). 

The grid size of GCMs is moving toward grids sufficiently fine to 
explicitly resolve many cloud systems, but the computational cost 
is enormous and, because of the importance of unresolved 
processes at still finer scales, convergence is by no means assured. A 
CRM can simulate clouds at meter- to kilometer-scale grid resolu- 
tions. Computational infrastructure typically limits the simulation 
of clouds and cloud systems by CRMs to a relatively small domain 
(<10 3 -km x 10 3 -km) and short time periods (<1 month). 
Grabowski (2001) and Khairoutdinov and Randall (2001) first 
proposed the use of 2D cloud-resolving models as a “super- 
parameterization" to simulate cloud processes within GCM grid- 
cells, replacing cumulus parameterizations. Arakawa (2004) 
describes this configuration as a multi-scale modeling framework 
(MMF). In the MMF, a non-hydrostatic 2D CRM takes the place of 
the single-column cumulus parameterization used in conventional 
GCMs (Randall et al., 2003a; Arakawa, 2004; Tao et al., 2009). 

There are two teams developing MMFs, a newer effort by God- 
dard and a longer running effort by Colorado State University (CSU). 
The CSU MMF combines the Community Atmosphere Model 3.0 
(CAM, Collins et al., 2006), the System for Atmospheric Modeling 
(SAM, Khairoutdinov and Randall, 2003), and the Community Land 
Model (CLM, Dai et al., 2003) to form the super-parameterized CAM 
(SP-CAM). The GCMs at the core of the Goddard and CSU MMFs 
share a common ancestor, the National Center for Atmospheric 
Research (NCAR) Community Climate Model Ver. 3 (CCM3), but 
underwent separate additional development by Goddard and NCAR 
researchers. 

Taking on phenomena that have been identified as difficult for 
GCMs to reproduce well, CSU researchers have shown better results 


with SP-CAM in reproducing the diurnal cycle of convection 
(DeMott et al., 2007), orogenic propagating cloud systems 
(Pritchard et al., 2011), subtropical low cloud fields (Blossey et al., 
2009), and precipitation anomalies associated with the 
Madden-Julian oscillation (Benedict and Randall, 2009) and El 
Nino— Southern Oscillation (ENSO, Khairoutdinov et al., 2008). Tao 
et al. (2009) compare the SP-CAM and an earlier version of the 
Goddard MMF. Both MMFs resulted in better representation of 
global energy and water cycles compared to GCMs with cumulus 
parameterizations but had their own set of biases from using 2D 
CRMs and prescribed sea surface temperatures. Researchers at 
Pacific Northwest National Laboratory (PNNL) added to SP-CAM an 
explicit-cloud parameterized-pollutant approach that links aerosol 
and chemical processes on the large-scale grid with statistics of 
cloud properties and processes resolved by the CRM (Wang et al., 
2011a, 2011b). The PNNL MMF can be used to study aerosol 
effects on cloud microphysics (indirect effect) globally, a study topic 
typically confined to a CRM-sized domain. 

In Tao et al. (2009), the differences between the CSU and God- 
dard MMFs were smaller than their differences with standard 
GCMs. These differences may become larger less from the diverging 
evolution of their parent GCMs than from the addition of additional 
model components. The emphasis here at Goddard on land/atmo- 
sphere interactions and hydrologic model development has 
produced an MMF in which “multi-scale” includes the land surface 
and a significant range of options are available to the user through 
LIS. Here, we describe the development and operation of the 
current Goddard MMF— LIS, focusing on the model coupling and its 
initial testing, particularly with respect to surface variables. 

3. Components of MMF— LIS 

Fig. 1 depicts the integration and coupling of the components of 
MMF— LIS. The three principal components are a finite-volume (fv) 
GCM (Goddard Earth Observing System Ver. 4, GEOS-4), a 2D CRM 
(Goddard Cumulus Ensemble, GCE), and the LIS, representing the 
large-scale atmospheric circulation, cloud processes, and land 
surface processes, respectively. All numerical analysis is written in 
FORTRAN90. The C language is used to expand object-oriented 
features already in FORTRAN90, providing a virtual object- 
oriented programming environment managing operations within 
and between components. The MMF— LIS components represent 
the work of several different teams of scientists and engineers at 
Goddard. The fvGCM was developed in the former NASA Data 
Assimilation Office. The successor to the Data Assimilation Office, 
the Global Modeling and Assimilation Office (GMAO), adopted the 
fvGCM as their first operational model. The CRM was developed in 
the Mesoscale Atmospheric Processes Laboratory, and the LIS was 
developed by the Hydrological Sciences Laboratory. Assisted by the 
Hydrological Sciences Laboratory, the Mesoscale Atmospheric 
Processes Laboratory performed the integration and coupling of the 
three MMF— LIS components. 

3.1. The Goddard Earth Observing System Ver. 4 (GEOS-4) 

The fvGCM of MMF— LIS, the GEOS-4, was constructed by 
combining the finite-volume dynamical core developed at Goddard 
(Lin, 2004) with the physics package of the NCAR CCM3, which 
represents a well-balanced set of processes with a long history of 
development and documentation (Kiehl et al., 1998). The unique 
features of the finite-volume dynamical core include an accurate 
conservative flux-form semi-Lagrangian transport algorithm with 
a monotonicity constraint on sub-grid distributions that is free of 
Gibbs oscillation (Lin and Rood, 1996, 1997), a terrain-following 
Lagrangian control-volume vertical coordinate (c-coordinate), 
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Fig. 1 . Integration and coupling of MMF-LIS. Black arrows depict the flow of forcing 
input from the fvGCM to GCE and LIS, and green arrows depict the flow of output from 
GCE and LIS to the fvGCM. Table 1 lists the specific input forcing variables to GCE and 
LIS. in the top and bottom left example maps, global daily pressure fields (hPa, black 
contours) and daily surface precipitation (mm, purple shading) are depicted. Global 
land cover is depicted in the bottom right map. 


a physically consistent integration of the pressure gradient force 
for a terrain-following coordinate (Lin and Rood, 1997; Lin, 1998), 
and a mass-, momentum-, and total-energy-conserving vertical 
remapping algorithm. 

Depending upon the application, the number of levels used in 
GEOS-4 varies between 32 and 64, while the horizontal grid 
spacing can vary between 2.5° and 0.125°. The GEOS-4 and its 
successor, the GEOS-5, have been applied in climate simulation, 
data assimilation, and weather prediction modes. Atlas et al. (2005, 
2007) and Shen et al. (2006a, b, 2010) tested their capability to 
simulate Atlantic hurricanes, adequately resolving problems like 
erratic track, abrupt re-curvature, intense extratropical transition, 
multiple landfall and re-intensification, and interaction among 
vortices. 

3.2. The Goddard Cumulus Ensemble (GCE) 

The GCE model has the longest history of the MMF— LIS 
components. It has been developed and refined at Goddard over 
two and a half decades for simulating convective clouds and cloud 
systems. The initial development and core features of the GCE 
model are in Tao and Simpson (1993), with reviews of the appli- 
cation of the GCE model to understanding precipitation processes 
in Tao (2003) and Tao et al. (2003). The GCE is a non-hydrostatic, 
anelastic CRM composed of prognostic equations for momentum, 
potential temperature, and water vapor mixing ratio. It includes 
solar and infrared radiative transfer processes, a Kessler-type 
two-category (cloud drops and rain) liquid water scheme, and 
a three-category (cloud ice, snow, and graupel/hail) bulk ice 


microphysics scheme. The GCE model’s default bulk microphysical 
scheme has been modified to improve graupel concentrations in 
the stratiform region and cloud ice concentrations in the upper 
troposphere (Lang et al., 2007, 2011 ). These changes better address 
saturation issues in these regions and result in more realistic 
column ice water contents for longer-term simulations (Zeng et al., 
2008, 2009). 

3.3. The Land Information System (LIS) 

The LIS is a scalable land data assimilation system that integrates 
a suite of advanced LSMs, high resolution satellite and observational 
data, data assimilation and parameter optimization techniques, and 
high-performance computing tools. The initial development and 
principal features of LIS are described in Kumar et al. (2006) and 
Peters-Lidard et al. (2007). The LIS infrastructure has three layers, 
1 ) a core that controls the program execution and I/O and manages 
the user-defined components; 2) a middle abstractions layer con- 
sisting of generic representations of LSMs, domains, land surface 
parameters, and running mode; 3) an extension of the abstractions 
layer for the user-selected LSM, domain, parameter set, and running 
mode. Executing LIS generates spatially and temporally distributed 
estimates of land surface processes using either observed or model- 
derived meteorology to constrain and force the user-specified 
LSMs. These models include the operational LSMs: CLM, Noah, and 
Catchment. In addition to new LSMs, multiple versions of the same 
LSM may be added and used within LIS. These features facilitate 
ensemble modeling studies and the benchmarking and sensitivity 
testing of new parameterizations and coupled modeling system 
configurations. The data assimilation functions of LIS are not 
currently used in MMF— LIS. Planned future upgrades will add data 
assimilation to the fvGCM, allowing future users to utilize the 
assimilation functions already in LIS. 

The three available running modes in LIS, analysis, forecast, and 
coupled, are described in Kumar et al. (2006, 2008). Initially, the 
coupled running mode was designed for local- to regional-scale 
modeling. Input land cover datasets masked out areas containing 
inland water bodies, wetlands, and glaciers, and parameterizations 
for evaporation from these surfaces were by-passed. During 
MMF— LIS development, these parameterizations were re-integrated, 
and a global input dataset with these land covers was created. 

4. Coupling and execution of MMF-LIS 

4.1. The atmosphere, fvGCM and GCE 

In Fig. 2, the atmosphere consists of an array of GCM gridcells 
each with a series of embedded 2D GCE models, representing 
a series of x-z slices of the atmosphere. The size of the GCM grid- 
cells, the number of embedded GCE models per gridcell, and the 
number of GCE internal gridcells are user-specified. Initialization, 
execution, and finalization of GCE operations are controlled by the 
fvGCM. An internal state variable (ISV) data structure is used to pass 
forcing and output variables such as temperature and humidity to 
and from the MMF-LIS components. The ISV has 5 degrees of 
freedom; 2 absolute spatial coordinates, latitude and longitude, and 
three relative spatial coordinates. For GCE, the latitude and longi- 
tude are the central geographical location of each fvGCM gridcell, 
and three Cartesian coordinates (x,y,z) designate the relative posi- 
tion of GCE internal gridcells. The ISV makes it possible to preserve 
both grid geographical and processor layout seamlessly. 

The passing of forcing ISVs from fvGCM to GCE initiates GCE 
execution. The atmospheric forcing variables passed to GCE such 
as temperature, humidity, and ozone are fvGCM-gridcell layer 
means, e.g., T(lon, lat, 0, 0, a). A coordinate transformation from the 
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Fig. 2. The black squares are an array of GCM gridcells and the red lines inside them 
are embedded 2D time-variant CRMs. For comparison, the black dots represent ID 
time-invariant cumulus parameterizations within a traditional GCM. Figure adapted 
from Tao et al. (2009). 


terrain-following fvGCM (7-coordinate system to the fixed GCE 
Cartesian coordinate system completes the transfer from the 
fvGCM to GCE. For land surface fluxes, there are two options, 1 ) an 
fvGCM-gridell mean that is then randomized across the GCE 
internal grid to approximate surface turbulence, 2) the perturba- 
tion occurring at each GCE internal gridcell, e.g., H'(lon, lat, x , 0, 0). 

During GCE execution, GCE mean layer temperature, humidity, 
and cloud properties such as cloud ice mixing ratio are calculated 
for each member of the GCE array from its internal perturbations. 
After execution of the GCE array, temporal averaging takes place to 
account for the fvGCM clock time elapsed during execution of the 
GCE array. The temperature and humidity tendencies for the fvGCM 
gridcell {e.g., [9T(lon. lat. 0, 0.(7) /3t] moist } are calculated from the 
array of GCE layer means, a Cartesian to rr-coordinate trans- 
formation, and the elapsed clock time. These results are then 
passed back to the fvGCM. There are two options for subsequently 
forcing LIS, 1 ) the fvGCM-gridell means at the lowest fvGCM model 
layer, 2) temporally-averaged GCE internal perturbations at the 
lowest GCE model layer. 

4.2. The land— atmosphere coupling with US 

The coupling of the fvGCM to LIS is enabled by the NASA Earth 
System Modeling Framework (ESMF, Hill et al„ 2004). The ESMF is 
open-source software to assist the development of high- 
performance, multi-component numerical earth system models. It 
offers a variety of data structures for transferring data between 
components, tools and utilities to ensure component interopera- 
bility with consistent component behavior, and libraries for re- 
gridding, time advancement, and other common modeling func- 
tions.The MMF— LIS coupling uses the ESMF Coupler Component 
(ESMF_CplComp), import and export state objects (ESMF_State), 
and the Time Manager utility (c./, Collins et al„ 2005; Balaji et al., 
2011). The Coupler is used only for transferring data between the 
fvGCM and LIS, as ESMF components and utilities are designed for 
data transfer and re-gridding, not physical/dynamical computations. 
Because the original design of LIS wrapped the ESMF Coupler 
superstructure around LIS (Kumar et al., 2006), the Coupler did not 
have to be added, only modified to recognize the ISVs used in the 
atmospheric components as forcing/output data container variables 
that could be pointed to by the import/export ESMF state objects. 

To initiate LIS execution, the Coupler takes in the import state 
object pointing to the ISVs of forcing data for LIS (Table 1 ) and maps 
the forcing data to the LIS grid. After LIS execution, any required 
temporal and spatial averaging takes place within LIS. An export 
state object then points to the output that will be mapped as fvGCM 


and GCE lower boundary conditions. The Time Manager accounts 
for the time elapsed during LIS operations. Because the MMF— LIS 
coupling is based on the template developed for coupling LIS to 
the mesoscale Weather Research and Forecasting (WRF) model, the 
sequence of interactions between fvGCM, the Coupler, and LIS is 
equivalent to the sequence for LIS— WRF diagrammed and 
described in Kumar et al. (2008). 

The first attempt to add two-way land surface processes to MMF 
involved coupling a single LSM, CLM 2.1, to fvGCM, in a horizontal 
structure, GCE <-» fvGCM <-» CLM. This was the MMF configuration 
used in Tao et al. (2009). Although CLM permits up to 7 surface 
types to be specified within its gridcells to approximate the 
underlying land surface heterogeneity, its limited I/O capabilities 
required spatial averaging such that both input forcing and output 
boundary conditions pass to and from fvGCM at fvGCM resolution, 
e.g., H(lon. lat, 0,0,0). The horizontal structure of the older MMF 
was adapted for MMF— LIS, swapping LIS and the ESMF Coupler for 
the single LSM, creating GCE <-> fvGCM <-> Coupler <-> LIS. 

A second strategy was considered, a vertical structure, fvGCM <-> 
GCE «-> Coupler «-» LIS, in which there would be multiple GCE-LIS 
couples embedded within each fvGCM gridcell. The vertical struc- 
ture would avoid the communication overhead of passing ISVs 
through fvGCM. However, it would complicate processor distribu- 
tion in parallel computing operations and incur significant compu- 
tational cost to start and stop many instances of LIS and organize the 
ISVs for the surface boundary conditions. In comparison, the hori- 
zontal strategy involves a single instance of LIS and a main LIS grid at 
the same resolution as fvGCM. For simulating higher resolution 
interaction between the atmosphere and the land surface, multiple 
surface types, “tiles”, at 1— lOMcm, can be specified within the LIS 
gridcells and the I/O passed to and from the Coupler at tile resolu- 
tion. The computational cost is several times greater to calculate and 
pass surface fluxes at cloud-scale (<10-km) resolutions. Because 
most land surface processes are parameterized as ID phenomena 
and LIS time steps are on the order of several minutes, the total 
computational cost of running one instance of LIS, even with cloud- 
scale tiles, is nominal compared to executing the GCE arrays. 

5. First results: a comparison of different land surface model 
physics 

5.1. Physics differences, CLM version 2.0 vs. 2A 

A significant functional advantage of coupling to LIS over a single 
LSM is the suite of abstractions available in LIS that expands the user’s 


Table 1 

List of forcing inputs to GCE (left column) and LIS (right column). The LIS forcing 
height is fixed at the lowest fvGCM model layer. The ozone and SST state variables 
are derived from NOAA weekly Reynolds Optimum Interpolation SST Analysis Ver. 2 
(Reynolds et al., 2002) and an ozone product merging the NASA Upper Atmosphere 
Research Satellite (UARS) ozone measurements (Ziemke et al., 1998) and the 
Atmospheric Model Intercomparison Project 2 (AMIP2) ozone dataset (Kanamitsu 
et al., 2002), respectively. 


GCE forcing 

LIS forcing 

Local solar time, t(lat, Ion) 

0 3 (lat, Ion, z), imposed 
SST(lat, Ion), imposed 
Temperature, T(lat, Ion, z) 
Humidity, q(lat, Ion, z) 

Winds, u(lat. Ion, z), v(lat, Ion, z) 
Advection (lat, Ion) of T and q, 
adv(T), adv(q) 

Sensible and latent heat flux, 
H(lat, Ion), LE(lat, Ion) 

Local solar time, t(lat, Ion) 

Winds, u, v(lat, Ion) 

Rain and snow rate, r(lat, Ion) 

Air temperature, T a j r (lat, Ion) 
Specific humidity, q(lat, Ion) 

Sea level pressure, P( lat, Ion) 
Shortwave radiation, SWj n (lat, Ion) 

Longwave radiation, LWi n (lat,lon) 
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capability to test codes and compare simulations using an ensemble 
modeling approach with different LSMs and/or versions of those 
models. The testing of new models or new model physics in LIS is 
accomplished by the addition or removal of components in the third 
layer (the extension of abstractions layer) where codes for specific 
LSMs reside. If new LSMs are added, communication between the 
new model and lower layers is established using the templates for I/O 
contained in the second layer. The LSMs added to LIS are no longer 
limited by their native I/O and tiling capability as these functions are 
determined and controlled by the top-level core layer. 

The original version of CLM in LIS was version 2.0. Added to LIS 
was the modified version of CLM, version 2.1, the same version of 
CLM directly coupled to MMF in the configuration of Tao et al. 
(2009). In CLM 2.1, there are several important changes to the 
model physics involving the computation of atmospheric forcing 
height, vegetation temperature, canopy interception of precipita- 
tion, and the drag coefficient between the underlying soil 
(or canopy surface) and the canopy air. All four of these state 
variables are used in determining the surface heat fluxes supplied 
to the atmospheric components. The underlying assumptions, 
the formulation of the equations for each variable, and the solution 
methods vary between versions. For example, the canopy inter- 
ception in CLM 2.0 is calculated assuming a uniform precipitation 
rate over the gridcell compared to an assumption of a logarithmic 
probability distribution function for precipitation rates in CLM 2.1. 

Because LIS had never before been run globally in coupled mode 
and new land covers and model physics were introduced, it was 
critical in the model development process to benchmark the 


MMF— LIS against the MMF configuration of Tao et al. (2009). We 
ran and compared two simulations of MMF— LIS, one using CLM 2.0 
and the other using CLM 2.1, to verify that MMF— LIS can produce 
physically realistic results and assess how changes to the physics of 
land surface processes would integrate over space and time to 
affect regional to global atmospheric circulations. 

5.2. Surface diagnostics 

The comparison of the different versions of CLM, original (2.0) 
versus modified (2.1), involved running MMF— LIS for two years, 
2007—2008, using each version at horizontal resolutions of 
2° x 2.5° for fvGCM and 4-km for GCE, with time steps of 30-min 
(fvGCM), 10-s (GCE), and 3-min (LIS). Vertical resolutions were 30 
layers (terrain-following) for fvGCM and 32 layers (fixed) for GCE. 
We specified the same simple surface heterogeneity as the MMF of 
Tao et al. (2009), 7 tiles per LIS gridcell and option 1 for land surface 
fluxes (fvGCM gridcell means). Because of the additional physical 
and computational complexity of cloud-scale heterogeneity, 
development of the MMF— LIS required that we first benchmark 
MMF— LIS against MMF using the same simple heterogeneity, 
leaving the testing of more complex representations of surface 
heterogeneity to future work. 

Fig. 3 summarizes the annual daily mean 2-m air temperatures 
for CLM 2.0 (original, top panels), CLM 2.1 (modified, middle 
panels), and their difference (original— modified, bottom panels) 
for 2007 (left panels) and 2008 (right panels). In 2007, most of the 
differences between the original and modified versions of CLM are 





Fig. 3. Maps of the annual daily mean and difference 2-m air temperatures in Kelvins for CLM 2.0 and CLM 2.1 for 2007 (a-c) and 2008 (d-f) at the resolution of the fvGCM 
(2° x 2.5°). The top panels (a, d) are the original, CLM 2.0, results. The middle panels (b, e) are the modified, CLM 2.1, results. The bottom panels (c, f) are difference maps of CLM 
2.0-CLM 2.1. The scales on the right column apply to the maps on the left column as well. 
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positive, implying CLM 2.0 tended to produce warmer air 
temperatures than CLM 2.1. Differences greater than 1 I< occur 
across much of Russia and Canada and in smaller areas of India, 
Australia, the Arabian Peninsula, and Antarctica. Noteworthy 
negative differences occur in Greenland, far-eastern Siberia, 
southern Brazil, and the North African coast. The trends in 2008 
amplify the major trends of 2007, with a larger area of positive 
differences in Eurasia and Canada to Alaska and expanded areas of 
negative differences in Greenland and North Africa to the Middle 
East. Differences in Antarctica between 30°E and 150°E change 
sign from 2007 to 2008 as the pool of coldest air (<220 K) in CLM 
2.1 in 2007 (panel b) is smaller in 2008 (panel e) versus CLM 2.0 
(panels a and d, respectively). The positive differences in the 
northern middle to high latitudes are in vegetated areas, and 
negative differences in Greenland, North Africa, Antarctica, and 
far-eastern Siberia are in bare or ice covered areas. 

Fig. 4 depicts the shortwave downward flux (SDF), typically the 
largest input to the surface energy budget. Because the first-order 
determinants of SDF are latitude, cloud fraction, and the optical 
thickness of clouds, the maps in Fig. 4 reflect zonal differences 
between the two versions in cloud cover and thus the radiation 
input to the surface energy budget. The areas of SDF less than 
120 W m~ 2 over the northern high latitudes are noticeably larger in 
CLM 2.1 (panels b and e) than in CLM 2.0 (panels a and d). The 
differences between the two models (bottom panels) are propor- 
tionally much larger north and south of 30°, over 30% in Central 
Siberia, but only 10—15% in Central America in 2007 (panel c). 


During the daytime over vegetated surfaces, net radiation is 
dominated by the shortwave downward term. 

The drag coefficient is important in determining how much 
SDF penetrates through the canopy and reaches the ground. 
Comparing the SDF in 2008 to the ground heat flux (Fig. 5a) 
indicates how much each version of CLM partitioned this energy 
into the fraction conducted into the surface versus the fraction 
available at the surface for evapotranspiration. The ground heat 
flux of the original CLM 2.0 tends to be greater in those areas with 
positive temperature differences in Fig. 3 and less in those areas 
with negative differences. Although on diurnal time scales, 
warmer surface temperatures due to enhanced ground heat flux 
feedback to the atmosphere through enhanced sensible heat flux. 
The enhanced sensible heat flux occurs only as long as it takes the 
surface layer air temperatures to adjust upward. The areas of 
positive and negative differences in Figs. 3 and 5 coincide such 
that the long-term mean sensible heat fluxes (not shown) in CLM 
2.0 and 2.1 are similar in these regions. Because of the heat storage 
capacity of soil, ground heat fluxes and thus surface temperature 
(2008, Fig. 5b), both instantaneous and long-term means, may 
reflect significant memory of radiation inputs. For LSMs like CLM 

2.0 that have a constant drag coefficient, Zeng et al. (2005) found 
that too much radiation was partitioned into the ground in areas 
with sparse canopies (e.g., Canadian and Siberian tundras) leading 
to warm surface temperature biases. The drag coefficient in CLM 

2.1 is allowed to vary with the friction velocity and thus the 
canopy thickness, ameliorating this bias. 
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Fig. 4. Maps of the annual daily mean and difference shortwave downward flux in W m 2 for CLM 2.0 and CLM 2.1 for 2007 (a-c) and 2008 (d-f). The resolution and layout of this 
figure are the same as Fig. 3. 
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Fig. 5. The 2008 annual daily mean difference maps for a) ground heat flux in W m 2 b) surface temperatures in Kelvins. For the sake of brevity, only maps for 2008 are shown. 


5.3. Effects on general circulation 

Fig. 4c depicts a large positive difference in SDF in Mexico in 2007. 
Fig. 6 contains time series of the gridcell centered at 26°N 102.5° W for 
the first month of the MMF— LIS simulations. Both simulations, CLM 
2.0 and CLM 2.1, start with identical initial conditions. Differences 
between the two versions begin to emerge on 2 Jan. By 3 Jan, the 
difference in the daily maximum SDF is more than 300 W itT 2 
because the total cloud amount time series are at opposite ends of the 
y-axis. Although the SDF and cloud amount variability track each 
other closely, surface temperature variability in both versions tracks 
the variability in SDF/cloud amount closely only for the first two 
weeks. After two weeks, there is less variability in any of the variables 
in Fig. 6 and a smaller lag between changes in SDF/cloud amounts and 
changes in surface temperatures in the original CLM 2.0. Differences 
between the drag coefficients and vegetation temperatures in the 
models contribute to the rapidly evolving differences between the 
two versions in the surface variables of Fig. 6. 

The gridcell examined in this section is an important source 
region for dry, warm air that, if it penetrates the Central US, forms 
an elevated mixed-layer over humid airmasses transported from 
the Gulf of Mexico. Together, these two air streams have long been 
recognized for their impact on the convective weather environ- 
ment in the Central US, particularly of the severe variety (Benjamin, 
1986; Lakhtakia and Warner, 1987; Lanicci et al., 1987). Trier et al. 
(2011) demonstrate how differences in model predictions of 
surface conditions in the high terrain of Mexico can result in 
significant differences in warm season precipitation in the Central 
US. In Fig. 7, the map panels show how quickly the simulated US 
regional circulations diverge after initialization. Subtle differences 
between the two simulations appear by Day 5 at 200 mb and Day 
10 for the sea level pressures. By Day 15, there are substantial 
differences in the locations, orientations, and magnitudes of all 
cyclones/troughs and anticyclones/ridges in the region. After Day 
16 in the critical Mexican gridcell (panels d and h), the sea level 
pressure time series of the original CLM 2.0 simulation indicates 
a series of transitory anticyclones versus a persistent longwave 
trough in the modified CLM 2.1 simulation. The anticyclones 
contribute to and are strengthened by the surface heating indicated 
in Fig. 6 by the CLM 2.0 simulation. The differences in the regional 
pressure patterns due to feedbacks between surface conditions, 
cloud amounts, and the upper troposphere take approximately 
10—14 days to emerge, consistent with the deterministic predict- 
ability limit of two weeks (reviews in Lewis, 2005; Yoden, 2007). 

6. Comparison to global gridded datasets 

6.1. Surface fluxes 

The FLUXNET is a global network of more than 500 microme- 
teorological tower sites (Baldocchi et al., 2001). The siting of these 


towers in a wide variety of biomes has made it possible to create 
gridded flux products useful for model validation from local to 
global scales (Jung et al., 2009, 2010; Blyth et al., 2010; Schlosser 
and Gao, 2010). In addition to comparing the MMF— LIS results to 





Fig. 6. Time series of the gridcell centered at 26°N 102.5°W (Mexico) for the first 
month of the MMF-LIS simulations for a) shortwave downward radiation flux, b) total 
cloud amount, c) surface temperature. Both simulations, original CLM 2.0 (black line) 
and modified CLM 2.1 (red line), start with identical initial conditions. 
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Fig. 7. North American regional maps of 200 mb temperatures (a-c) and sea level pressures (d-g) for the CLM 2.0 (black line) and CLM 2.1 (red line) simulations after 5 days (a, e), 
10 days (b, f), and 15 days (c, g). The bottom panels are the time series for the gridcell centered at 26°N 102.5°W (Mexico) for d) 200 mb temperature and h) sea level pressure. 
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a FLUXNET gridded product, we chose the NASA GMAO reanalysis 
product, the Modern Era Retrospective-analysis for Research and 
Applications (MERRA). The MERRA is based on the Goddard Earth 
Observing System Ver. 5 (GEOS-5, Bosilovich et al„ 2006, 2011 ). The 
GEOS-5 has the same dynamical core as the GEOS-4 but improved 
moist process parameterizations and data assimilation (Rienecker 
et al., 2007). The data assimilation system of GMAO’s interactive 
land-ocean-atmosphere modeling framework is NCEP’s Gridpoint 
Statistical Interpolation (Wu et al., 2002). By comparison, the 
MMF— LIS results are from free-running simulations, i.e„ lacking 
a data assimilation system, and partially interactive, using observed 
SST/sea ice rather than a full-physics ocean model. 

Fig. 8 depicts June— August 2007 mean latent heat flux for 
FLUXNET, MERRA, and the two MMF— LIS simulations. Compared 
to the FLUXNET and to MMF— LIS, the MERRA map has significantly 
more latent heat flux over the humid regions of Central Africa, 
Central to South America, East Asia, Indonesia, and eastern North 
America. Reichle et al. (2011 ) attribute this positive bias in areas to 
excessive interception by dense canopies. The MMF— LIS maps 
have fewer and smaller areas with latent heat fluxes greater than 
120 W m -2 , although there is a notable hot spot in the extensive 
wetlands around the Ob River in Central Siberia that is not re- 
flected in the FLUXNET map. In arid Australia and Central Asia, the 
MERRA performs better than both MMF— LIS simulations. 
Comparing the different versions of CLM in MMF— LIS, the areas of 
elevated latent heat flux in Eurasia and North America are smaller 
and weaker in the modified CLM 2.1. The one area of latent heat 
flux greater than 120 W itT 2 in the FLUXNET map is in the 
Central US. It is better represented in the original CLM 2.0 map, 
although the modified CLM 2.1 performs better than the CLM 2.0 
in the important dry airmass source region in the Mexican high 


terrain and in the northern tundras affected by the warm bias in 
CLM 2.0. 

6.2. Precipitation 

There are a number of global gridded precipitation datasets 
available to assess model performance from diurnal time scales 
onward. The CMORPH (Climate Prediction Center MORPH ing 
technique) is composed of precipitation estimates derived from low 
Earth orbiting satellite microwave observations whose features are 
transported via spatial propagation information obtained from 
geostationary satellite IR data (J°y ce et al., 2004). Fig. 9 contains 
map views and latitudinal and longitudinal cross sections of 
CMORPH, MERRA, and the two MMF— LIS simulations of the mean 
daily precipitation for June— August 2007 (50°N-50°S). There are 
areas where all of the models produce too little precipitation (the 
equatorial Atlantic, northern Europe, Central US), but there is 
generally more precipitation in the model output than in the 
CMORPH. In the Indian-Asian monsoon region, all of the models 
produce a much larger area of heavy (>10 mm day -1 ) precipitation 
than in CMORPH. This is particularly acute over the Western Pacific 
warm pool in the MMF— LIS simulations versus the Indian monsoon 
region in MERRA. Around Indonesia and in northern South Amer- 
ica, Central Asia, and southeastern Australia MERRA is both close to 
the CMORPH and clearly better than both MMF— LIS simulations. 
Both MMF— LIS simulations outperform MERRA in sub-Saharan 
Africa. Over the Pacific, the modified CLM 2.1 outperforms 
MERRA east and north of the western warm pool. In the latitudinal 
cross section, the boreal summer rainfall peak is wider in the 
MERRA, contributing to a global mean rainfall rate 30% greater than 
CMORPH versus 23% greater for MMF— LIS. 




Fig. 8. Global maps of June- August 2007 mean latent heat flux in W m 2 for a) FLUXNET, b) CLM 2.0, c) MERRA, d) CLM 2.1. The color scale is the same for all maps. Grid resolution 
of each map is indicated. 
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Fig. 9. Global mean daily precipitation in mm day 1 for June-August 2007 (50°N-50°S) with map views and latitudinal and longitudinal cross sections of a) CMORPH, b) CLM 2.0, c) 


MERRA, d) CLM 2.1. 


There are important differences between the MMF— LIS simu- 
lations. The modified CLM 2.1 compares better to CMORPH than the 
original CLM 2.0 in most of the Pacific, the South Atlantic, northeast 
Asia, the eastern and Central US, and sub-Saharan Africa. The areas 
affected by the Mexican high terrain source region, the Central US 
and Mexico, are wetter in the modified CLM 2.1 map. The original 
CLM 2.0 performs better than the modified CLM 2.1 in a couple 
much smaller areas, the northwestern Atlantic along the US coast, 
where CLM 2.1 is too dry, and in southeastern Australia, where CLM 
2.1 is too wet. In the South Atlantic, the modified CLM 2.1 has 
a lower precipitation maximum than MERRA, but it is higher than 
the original CLM 2.0, and in none of the models is it tangential to 
the South American coast as in CMORPH. Comparing the longitu- 
dinal cross sections, the peak at 150°E is lower in the modified CLM 
2.1, contributing to the 1% reduction in global mean rainfall versus 
the original CLM 2.0. The results in Fig. 9 for the modified CLM 2.1 
are highly similar to the results from the MMF configuration of Tao 
et al. (2009), completing a successful benchmarking of the 
MMF— LIS code against the older configuration. 

In a study of precipitation biases over the tropical Atlantic by 6 
different free-running GCMs by Biasutti et al. (2006), the cumulus 
parameterizations in the models overestimated the correlation 
between convective precipitation and surface humidity and 
underestimated the correlation to upper tropospheric humidity. 
They were thus acutely sensitive to SST distribution and latent heat 
flux. This sensitivity is compounded in both parameterized and 
resolved convection schemes by several additional factors. These 
schemes produce unrealistically high updraft and downdraft mass 
fluxes and precipitation efficiencies, raining at column humidities 


lower than observed, particularly in humid regions (Sui et al., 2007; 
Kain et al., 2008; Weisman et al., 2008). Unable to sustain high 
column humidity, they may rain prematurely, affecting the simu- 
lation of the diurnal cycle (Mohr et al., 2003; Tompkins and Di 
Giuseppe, 2003; Bernie et al., 2007) and of convective develop- 
ment associated with fronts and troughs (Zeng et al., 2007; 
Weisman et al., 2008). The timing and amount of precipitation in 
simulations with resolved convection is also highly sensitive to the 
choice of ice microphysical schemes (Li et al., 2009a, b; Satoh et al., 
2010) and to the use of cyclic lateral boundary conditions that trap 
convection within the CRM, producing artificially long lifetimes 
(Tao et al., 2009). 

For coupled ocean-atmosphere model configurations, adjust- 
ment between the surface atmospheric and oceanic layers tends to 
warm surface air temperatures and reduce total cloud amount 
(Costa et al., 2001; Biasutti et al., 2005; Raisanen et al., 2008). This 
effect occurs even in slab ocean models, although anomalous latent 
heat fluxes may develop in slab models to compensate for imposed 
ocean heat transport, an effect not present in full-physics mixed- 
layer ocean models (Sutton and Mathieu, 2002). Although the use 
of mixed-layer ocean models and data assimilation systems can 
reduce precipitation biases, the results in Figs. 8 and 9 and of 
previous studies (e.g„ Biasutti etal., 2006; Shutts and Palmer, 2007; 
Stan et al., 2010) suggest improved model physics is still important. 
Issues with model macro- and microphysics, surface observational 
networks, satellite radiance retrieval schemes, grid resolution, and 
the expenses associated with computing and storage remain major 
issues in global modeling (reviews in Zhang et al., 2005; Guilyardi 
et al., 2009). 
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7. Conclusions 

7.1. Summary 

The three principal components of the MMF— LIS are a finite- 
volume GCM (Goddard Earth Observing System Ver. 4, GEOS-4), 
a 2D CRM (Goddard Cumulus Ensemble, GCE), and the Land Infor- 
mation System (LIS), representing the large-scale atmospheric 
circulation, cloud processes, and land surface processes, respec- 
tively. This global multi-scale modeling framework exchanges the 
single-column cumulus parameterization common to conventional 
GCMs for a 2D CRM running at non-hydrostatic resolutions 
(<10-km). The MMF— LIS has a horizontal coupling strategy, GCE <-> 
fvGCM <-> Coupler <-» LIS, that uses an Earth System Modeling 
Framework Coupler Component to transfer I/O between fvGCM and 
LIS. This strategy offers significant computational efficiency by 
simplifying processor distribution during parallel operations and 
minimizing I/O overhead while maximizing its flexibility. The I/O 
between the atmospheric and land surface components can be 
passed either as fvGCM gridcell means or as cloud-scale (1— 10-km) 
tiles. The latter option permits cloud-scale two-way interaction at 
a nominal additional computational cost. Coupling with LIS versus 
a single LSM expands the user’s ability to perform ensemble 
modeling with different LSMs and add and test new model physics. 

We tested two different versions of CLM, 2.0 and 2.1, simulating 
2007—2008 to assess whether the new MMF— LIS could produce 
physically realistic results and assess how changes in model physics 
might affect regional atmospheric dynamics. Our benchmark was 
the MMF configuration of Tao et al. (2009) that ran CLM 2.1 with 
simple surface heterogeneity (7 tiles per LIS gridcell). Comparing 
mean annual global surface temperature maps, large areas of the 
northern hemisphere were several degrees warmer in the CLM 2.0 
simulation. Positive feedbacks between warmer surface, warmer 
soil, and lower cloud cover set and reinforced the warming 
tendency, making the positive (original— modified) differences 
greater in area and magnitude in 2008. By starting each simulation 
with the same initial conditions, we examined how quickly differ- 
ences between both simulations would emerge and affect the 
regional North American circulation. Differences emerged after 2 
days for surface temperatures, 5 days for upper tropospheric 
temperatures, and 10 days for the synoptic-scale surface pressure 
fields. By the end of the deterministic prediction limit of two weeks, 
all of the synoptic-scale cyclones and anticyclones in North America 
had different positions, orientations, and magnitudes. 

The first results from MMF— LIS benchmarked successfully 
against an older configuration of MMF. Examining two MMF— LIS 
simulations with different LSM physics revealed sensitivity to 
changes in model physics, producing noteworthy differences in 
both surface variables and regional circulations. From free-running 
simulations, the surface fluxes and precipitation over land from 
MMF— LIS compared well to global gridded observational and 
reanalysis (MERRA) datasets. The spatial coverage of precipitation 
greater than 2 mm day -1 from MMF— LIS was comparable to 
observed. The weakest aspect of the MMF— LIS simulations of global 
precipitation was the positive precipitation bias in the Indian-Asian 
monsoon region. Comparing the CLM 2.1 precipitation and latent 
heat flux maps to the CLM 2.0 simulation and observations showed 
improvement in many areas, particularly in the northern hemi- 
sphere where CLM 2.0 had a notable warm bias. 

7.2. Future directions in global land/atmosphere modeling 

Future improvements in representing land/atmosphere inter- 
actions in operational global forecasting models includes 
improving grid resolutions, better individual physical process 


representations, and data assimilation, both from better assimila- 
tion methodologies and additional data streams (Saha et al., 2010; 
Bosilovich et al., 2011; Riddaway, 2011; Werth and Garrett, 2011). 
Adding multi-LSM capability through LIS or a LIS-like interface is 
not mentioned. Reducing GCM grid resolution will improve the 
representation of surface heterogeneity by underlying LSMs and 
thus land/atmosphere interactions. Below 10-km, the hydrostatic 
approximation breaks down, requiring new model formulations in 
operational agencies that accommodate the potentially opposing 
requirements of computational speed and physical accuracy while 
limiting code development time and effort where possible. 

At the ECMWF, work is on-going to develop a non-hydrostatic 
global model that, as computational infrastructure allows, 
permits the explicit resolution of moist convection with land/ 
atmosphere interactions at a comparable scale (Yessad and Wedi, 
2011). To save development effort, the existing dynamical core of 
the Meteo-France’s mesoscale forecasting model Aire Limitee 
Adaptation Dynamique Developpement International (ALADIN) is 
being evaluated for this role (Benard et al., 2010). At NCEP, separate 
operational execution of WRF (with Noah) at 4-km with explicit 
convection takes place but only over the continental United States 
and its outlying states and territories (Kain et al., 2008; Carley et al., 
2011). Although this area will not be expanded in the near term, 
a single nested structure is being created from the physical and 
dynamical cores of NCEP’s large-scale Eta and mesoscale WRF 
models linked by ESMF components and interfaces (Iredell and 
Black, 2011). 

The Goddard MMF— LIS is currently upgrading to GEOS-5. The 
GEOS-5 is significantly computationally more efficient and will add 
new functionality to MMF— LIS, in particular, a non-hydrostatic 
running option, a fully coupled deep-ocean component, and data 
assimilation capabilities that will allow future utilization of land 
and atmosphere data assimilation. The latter options will address 
the positive precipitation biases over the tropical oceans that 
feedback to the land/atmosphere interaction. The next step 
conceptually for the present generation of MMF configurations that 
use zonal arrays of CRMs is described by Arakawa et al. (2011) as 
a “Quasi-3D MMF” using a network of CRMs forming a grid with 
large gaps. The gap sizes would be determined by available 
computational infrastructure and desired application. As compu- 
tational infrastructure evolves and gaps shrink, a true 3D MMF will 
emerge. 
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